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Abstract 

The  gas  diffusion  layer  of  a  polymer  electrolyte  membrane  (PEM)  fuel  cell  is  a  porous  medium  generally  made  of  carbon  cloth  or  paper.  The  gas 
diffusion  layer  has  been  modeled  conventionally  as  a  homogeneous  porous  medium  with  a  constant  permeability  in  the  literature  of  PEM  fuel  cell. 
However,  in  fact,  the  permeability  of  such  fibrous  porous  medium  is  strongly  affected  by  the  fiber  orientation  having  non-isotropic  permeability. 
In  this  work,  the  lattice  Boltzmann  (LB)  method  is  applied  to  the  multi-phase  flow  phenomenon  in  the  inhomogeneous  gas  diffusion  layer  of  a 
PEM  fuel  cell.  The  inhomogeneous  porous  structure  of  the  carbon  cloth  and  carbon  paper  has  been  modeled  as  void  space  and  porous  area  using 
Stokes/Brinkman  formulation  and  void  space  and  impermeable  fiber  distributions  obtained  from  various  microscopic  images.  The  permeability 
of  the  porous  medium  is  calculated  and  compared  to  the  experimental  measurements  in  literature  showing  a  good  agreement.  Simulation  results 
for  various  fiber  distributions  indicate  that  the  permeability  of  the  medium  is  strongly  influenced  by  the  effect  of  fiber  orientation.  Present  lattice 
Boltzmann  flow  models  are  applied  to  the  multi-phase  flow  simulations  by  incorporating  multi-component  LB  model  with  inter-particle  interaction 
forces.  The  model  successfully  simulates  the  complicated  unsteady  behaviors  of  liquid  droplet  motion  in  the  porous  medium  providing  a  useful 
tool  to  investigate  the  mechanism  of  liquid  water  accumulation/removal  in  a  gas  diffusion  layer  of  a  PEM  fuel  cell. 

©  2007  Elsevier  B.Y.  All  rights  reserved. 
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1.  Introduction 

Owing  to  its  high  efficiency,  zero  emission  and  low  operating 
temperature,  polymer  electrolyte  membrane  (PEM)  fuel  cell  has 
been  widely  considered  as  a  promising  alternative  power  source 
for  mobile,  portable  and  stationary  co-generation.  A  PEM  fuel 
cell  consists  of  a  membrane  electrolyte  assembly  (MEA)  sand¬ 
wiched  between  two  bipolar  plates  as  shown  in  Fig.  1.  The  gas 
diffusion  layer  (GDL),  catalyst  layer  and  polymer  electrolyte 
membrane  are  referred  to  as  MEA  where  current  is  produced. 
The  GDL  is  a  porous  medium  generally  made  of  carbon  cloth  or 
paper  and  provides  essential  roles  in  a  PEM  fuel  cell  operation 
such  as  electronic  conductivity,  reactant  access  to  the  catalyst 
layers  and  product  removal  from  it.  The  performance  of  a  PEM 
fuel  cell  deteriorates  due  to  liquid  water  accumulated  in  the  GDL 
and,  therefore,  optimizing  the  key  parameters  of  the  GDL  such  as 
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thickness,  porosity,  permeability  and  wetting  characteristics  has 
been  of  considerable  importance  for  the  PEM  fuel  cell  design, 
operation  and  eventual  commercialization. 

Studies  on  the  GDL  have  been  led  mostly  by  experimental 
studies  [1-5].  The  pores  can  be  hydrophobic  with  optimal  value 
of  polytetrafluoroethylene  (PTFE)  not  to  be  congested  with  liq¬ 
uid  water  [1].  Water  flow  in  GDL  is  experimentally  investigated 
for  various  thicknesses  and  pore  sizes  in  Benziger  et  al.  [2]. 
Paganin  et  al.  [1]  observed  a  performance  decrement  at  higher 
current  densities  when  the  GDL  thickness  is  increased.  The 
thickness  of  the  GDL  was  optimized  by  a  mathematical  model¬ 
ing  in  Inoue  et  al.  [3]  and  by  cell  performance  tests  in  Lee  et  al. 
[4]  and  Jordan  et  al.  [5],  with  various  GDL  parameters  such  as 
porosity  and  thickness.  A  thin  GDL  with  small  porosity  results 
in  good  electrical  conductivity,  however  efficient  mass  transport 
requires  large  pores,  i.e.  large  porosity.  With  a  small  and  long 
serpentine  flow  channel  that  can  be  considered  as  many  par¬ 
allel  channels  connected  in  series,  reactants  can  directly  cross 
to  neighboring  flow  channels  due  to  high  pressure  gradient  and 
short  flow  path  [6-9].  The  amount  of  cross  flow  is  strongly  cor- 
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Nomenclature 

Ca  capillary  number 

c  8x/8t 

cs  sound  speed 

e;  velocity  vector  for  the  lattice  Boltzmann  model 

fi(x,  t)  distribution  function 

/fq(x,  t)  equilibrium  distribution  function 
G  a  potential  that  couples  nearest  and  next-nearest 
neighbors. 

k  species  number 

K  permeability  of  the  medium 

Ktow  tow  permeability 

Ktow  permeability  tensor 

L  characteristic  length 

m  molecular  weight 

N  number  of  grid  points 

P  pressure 

R  universal  gas  constant 

ReL  Reynolds  number 

T  absolute  temperature 

u  velocity  vector 

ux  velocity  in  x-direction 

uy  velocity  in  y-direction 

V  characteristic  velocity 

X  position  vector 

Greek  letters 

ft  a  parameter  controlling  the  magnitude  of  the 

momentum  sink 

8X  the  distance  between  abutting  lattices 

8t  time  step 

/x  viscosity  of  fluid 

v  kinematic  viscosity 

p  density  of  fluid 

a  species  number 

r  relaxation  time 

0  porosity,  Nvoid/AWl 

t  an  arbitrary  function  of  density  distribution  for 

species  o 

cot  weight  coefficient  for  the  lattice  Boltzmann  model 

Subscripts 
in  inlet 

max  maximum  value 

min  minimum  value 

out  outlet 

void  void  region  in  the  GDL 

total  the  total  region  in  the  GDL 

x  x-direction 

y  y-direction 

z  z-direction 


related  with  the  thickness  and  permeability  of  the  GDL  in  a 
PEM  fuel  cell  with  a  serpentine  flow  channel  [6-8].  It  was  also 
shown  that  such  cross  flow  can  be  effective  for  the  removal  of 
liquid  water  from  the  gas  diffusion  layer  [9].  The  GDL  made 
of  carbon  cloth  or  paper  normally  has  anisotropic  permeabil¬ 
ity  for  two  major  orientations:  in-plane  and  through  the  plane 
direction.  The  in-plane  permeabilities  of  various  GDLs  were 
measured  from  the  order  of  10“ 12  to  10“ 11  m2  in  Feser  et  al. 
[10]. 

A  conventional  numerical  treatment  of  the  GDL  in  a  PEM 
fuel  cell  is  modeling  the  entire  GDL  as  a  homogeneous  porous 
medium  having  isotropic  or  anisotropic  permeability  [11-13]. 
The  method  is  numerically  efficient  since  it  does  not  demand 
a  fine  grid  to  capture  the  microscopic  structure  of  the  GDL. 
However  the  result  of  this  model  is  mainly  dependant  on  the 
given  value  of  permeability  which  may  not  be  known  for  the 
various  porous  materials  of  GDLs  and  it  is  not  proper  to  utilize 
the  method  to  investigate  the  effect  of  fiber  structure  and  sur¬ 
face  condition  on  the  permeability  of  the  medium.  Although  it 
may  require  large  computing  power,  the  true  geometry  of  the 
porous  medium  can  be  modeled  directly  from  the  microscopic 
photography  of  the  medium.  Joshi  et  al.  [14]  applied  the  lat¬ 
tice  Boltzmann  (LB)  method  to  a  modeling  of  multi-component 
gas  transport  in  a  solid  oxide  fuel  cell  (SOFC).  The  porous 
structure  of  the  SOFC  is  modeled  as  void  space  and  imper¬ 
meable  solid  obstacle  from  the  scanning  electron  microscope 
(SEM)  image  of  the  porous  electrode.  The  LB  model  could  suc¬ 
cessfully  simulate  the  diffusive  flow  in  the  porous  electrode  of 
SOFC  without  empirical  modification  of  diffusion  coefficients 
using  medium  porosity  and  tortuosity  [14].  An  intermediate 
scheme  to  model  the  flow  in  a  porous  electrode  may  be  the 
application  of  the  Stokes/Brinkman  formulations  to  allow  fluid 
transport  through  finite-sized  porous  objects  by  treating  them  as 
an  effective  medium  with  a  known  permeability  [15,16].  Spaid 
and  Phelan  [15]  showed  that  the  effective  permeability  of  the 
porous  medium  calculated  by  LB  model  well  matches  an  ana¬ 
lytical  calculation  by  Phelan  and  Wise  [17].  Park  et  al.  [16] 
performed  3D  LB  simulation  and  observed  that  the  permeabil¬ 
ity  of  porous  medium  is  strongly  dependant  on  the  fiber  tow 
orientation. 

The  LB  method  has  been  accepted  as  a  new  computational 
tool  for  a  variety  of  fluid  transport  phenomena  [18].  It  was 
applied  to  incompressible  fluid  flows  [19,20],  transport  of  pas¬ 
sive  scalars  [21],  miscible  and  immiscible  fluids  in  complex 
geometries  [22]  and  two-phase  flow  with  phase  change  [23] .  The 
kinetic  nature  of  the  LB  method  was  also  shown  to  be  applica¬ 
ble  to  simulation  of  chemical  reaction  in  micro-  and  meso-scopic 
flow  [24]  and  electrokinetic  transport  phenomena  [25].  In  this 
work,  we  have  modeled  the  two  representative  GDL  materials, 
woven  carbon  cloth  and  carbon  paper,  via  LB  method.  The  true 
porous  structure  of  carbon  paper  was  directly  taken  from  micro¬ 
scopic  images  following  Joshi  et  al.  [14].  The  permeability  has 
been  obtained  for  the  various  sample  images  of  carbon  papers 
and  compared  with  the  experimental  measurements  by  Feser  et 
al.  [10].  We  also  present  an  unsteady  two-phase  LB  flow  simu¬ 
lation  which  is  applied  to  the  study  of  liquid  water  removal  from 
the  porous  electrode  of  a  PEM  fuel  cell. 
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Membrane  Electrode  Assembly  (MEA) 


Fig.  1.  Schematic  of  a  polymer  electrolyte  membrane  (PEM)  fuel  cell  assembly. 


2.  Numerical  section 

2.1.  Lattice  Boltzmann  model 

The  LB  method  is  based  on  a  finite  number  of  identical  par¬ 
ticles  that  go  through  collision  and  propagation  successively  on 
prefixed  paths  in  space.  The  following  single-component  lat¬ 
tice  Boltzmann  equation  describes  evolution  of  the  distribution 
function, //(x,  t),  with  the  BGK  collision  term  [26,27]: 


where  c  =  8x/8t ,  and  8X  is  the  distance  between  lattice  points.  The 
speed  of  the  sound  is  c/\[3  in  the  phase  space  and  accordingly 
we  have  RT  =  c2/3  and  the  fluid  pressure  is  given  in  terms  of  the 
fluid  density  as  p  =  p/3.  The  macroscopic  number  density,  p(x, 
0,  and  the  velocity,  u(x,  t ),  of  the  fluid  are  obtained  as 

8 

p  =  ffmfi ,  (4) 

;=0 


fi(x  +  eiSt ,  t  +  St)~  fi(x,  t)  =  - 


fi(x,  Q  -  fCHx,  t) 

Ty 


8 

(1)  pu  =  ffmftet, 
i=  0 


(5) 


where  e/’s  are  the  discrete  velocities,  5/  is  the  time  step  and  zv  is 
the  relaxation  time.  /zeq  represents  the  equilibrium  distribution 
of  fi  given  as 


where  the  molecular  weight,  m,  is  assumed  to  be  unity  in  the 
present  work.  The  kinematic  viscosity  is  related  to  the  relaxation 
time  as 


fi\p,  M)  =  U>iP 


1  + 


RT 


u  (e;  ■  ur 

+ 


2  (RT)1 


u  •  u 
2RT 


(2rv  -  1)  8j 

6  8,' 


(6) 


4/9,  i  =  0 
1/9,  i=  1,2,  3,4 
1/36,  i  =  5,6,7,  8 


(2) 


where  oof  s  are  the  associated  weight  coefficients,  R  is  the  uni¬ 
versal  gas  constant  and  T  is  the  absolute  temperature.  In  Fig.  2 
the  velocity  vectors,  ez,  for  the  two-dimensional  9- speed  model 
(D2Q9)  are  shown  to  be 


(0, 0), 


(  i-1  i- 1  \ 

COS - 71,  Sin - 71  , 

V  2  2  J 


r-  (  i-  4—1/2  ./ -4-1/2  \ 

V2 c  cos - tv,  sin - tv  , 

V  2  2  ) 


for  i  =  0 
for  i  =  1,2,  3,  4 

for  i  =  5,  6,  7,  8 

(3) 


Through  Chapman-Enskog  expansion,  Eqs.  (1)  and  (2)  lead 
to  the  continuity  and  the  Navier-Stokes  equations  near  the 


(7)  (4)  (8) 


Fig.  2.  Lattice  structure  of  the  2D  (D2Q9)  lattice  Boltzmann  model. 
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incompressible  limit  [28]: 


2.3.  Boundary  conditions 


dp 

Yt 


+  V  •  (pu)  =  0, 


P 


—  +  (w.V)w 
ot 


—Wp  +  V  •  [pv(Vu  +  Vmt)]. 


(7) 

(8) 


2.2.  Multi-component  lattice  Boltzmann  model 

The  multi-component  LB  model  with  inter-particle  interac¬ 
tion  [20,29]  was  developed  for  simulation  of  multi-phase  flows 
and  phase  transitions.  When  the  interaction  is  weak  this  model 
can  be  applied  to  simulate  diffusion  due  to  various  driving  mech¬ 
anisms  [20,25,29].  If  any  external  force  is  applied  to  component 
k  the  momentum  should  be  incremented  correspondingly.  In  LB 
method,  this  is  simply  achieved  by  replacing  u  in  each  den¬ 
sity  equilibrium  function  with  u'  +  T0~¥0lp0.  General  form  of 
interaction  force  between  different  components  at  site  x  is 

Fsa  =  x  +  ei)fo(x  +  el)ei,  (9) 

i  o 

where  a  is  an  arbitrary  function  of  density  distribution,  taken 
as  fia  in  this  study.  In  the  D2Q9  lattice  model,  G0a  describes 
a  potential  that  couples  nearest  and  next-nearest  neighbors  as 
follows: 


The  practical  flow  situation  in  a  PEM  fuel  cell  employs  some 
forms  of  flow  regulations  at  the  inlet,  outlet  and  other  boundaries. 
Common  types  of  regulations  are  constant-velocity,  constant 
pressure  and  no-slip  wall.  For  the  D2Q9  model,  we  utilized  the 
bounce-back  rule  to  determine  unknown  part  of  particle  distri¬ 
bution  functions  [31].  As  an  example,  it  is  considered  that  the 
boundary  is  aligned  with  the  v-direction  with/4,/7,/8  pointing 
into  the  wall  as  shown  in  Fig.  2.  After  streaming, /o,/i,/3,/4, 
/7,/8  are  known.  Supposing  that  ux  and  uy  are  specified  on  the 
boundary  Eqs.  (4)  and  (5)  are  used  to  determine/2,/5,/6  and  p 
as  follows 

h  +  fs  +  h  =  P  -  (/o  +  /1  +  h  +  h  +  fi  +  h\  (15) 

fs  ~  h  =  pux  ~  (/1  -  h  ~  fi  +  h\  (16) 

h  +  fs  +  fa  =  Puy  +  (/4  +  fl  +  f%),  (17) 

P  =  -7^—Vfo  +  fl  +  h  +  2(/4  +  fi  +  m  (18) 

1  ~  Uy 

To  close  the  system  of  equations  above,  we  assume  the 
bounce-back  rule  for  the  non-equilibrium  part  of  the  particle 
distribution  normal  to  the  boundary  (in  this  case,  = 

/4  —  f^).  With/2  known  and/6  can  be  found,  thus 


Gc Ta(x,  Xf) 


ga&,  |*-*'|  =  1 

goo/ 4,  |X-X'|  =  V2  . 

0,  otherwise 


(10) 


Here  gGa  is  the  strength  of  the  inter-particle  potential  between 
component  a  and  a.  The  wall  is  regarded  as  a  phase  with  a 
constant  number  density.  The  interactive  force  between  the  fluid 
and  wall  is  described  as 


F wo  —  VG(x)^>w(*,  ^  _|_  ^)/w(x  +  €{)€i , 

i 


(ID 


where  the  /w  is  the  number  density  of  the  wall,  which  is  con¬ 
stant  at  the  wall  and  zero  elsewhere.  Ga w  is  similar  with  G0a 
while  positive  for  a  nonwetting  and  negative  for  wetting  fluid. 
Note  that  introduction  of  fluid-solid  interaction  has  no  effect  on 
the  macroscopic  equations  since  it  only  exists  at  the  solid/fluid 
interface.  The  Chapman-Enskog  expansion  procedure  can  be 
carried  out  to  obtain  the  following  continuity  and  momentum 
equations  for  the  fluid  mixture  as  a  single  fluid. 

%  +  V  •  (pu)  =  0  (12) 

ot 


p 


— —  +  (u  •  V)w 
ot 


— V p  +  V  •  [pv(Vu  +  Vmt)], 


(13) 


where  p  =  ^2aPo  is  the  total  density  of  the  fluid  mixture,  and 
the  whole  fluid  velocity  u  is  defined  as 


pu  =  Y\  PoUcj  +  -Fa. 
z — Jo  2 


(14) 


fl  =  U  +  ^PUy,  (19) 

fs  =  fl  -  i(/i  -  h)  +  +  ipw-y,  (20) 

h  =  h  -  i(/l  -  fi)  -  1 PUx  +  (21) 

The  collision  is  also  applied  to  the  boundary  nodes.  No  slip 
wall  can  be  achieved  applying  velocities  on  the  wall  equals  zero 
(ux  -uy-  0)  on  the  boundary.  Similarly,  when  the  density  (pres¬ 
sure)  is  known  on  the  boundary  and  the  velocity  is  assumed  to 
be  normal  to  the  boundary  ( ux  =  0),  the  unknowns  (uy,fs,fe)  can 
be  determined  from  Eqs.  (17)-(20)  as 


-  fo  +  fi  +  h  +  2(A  +  fi  +  h) 

Uy  =  1 - 

Pin 

2 

fi  =  U  + 

fs  =  fi  ~  i(/i  -  fi)  + 
h  =  h~  i(/i  -  fi)  +  ^puy. 


(22) 

(23) 

(24) 

(25) 


The  entire  operation  is  executed  after  streaming  and  before 
collision.  Obstacles  (carbon  fiber)  within  the  numerical  domain 
are  also  treated  like  impermeable  walls  with  no  mass  flux  in  a 
direction  normal  to  their  surface.  The  tangential  velocity  of  all 
species  at  surface  of  the  solid  obstacles  is  assumed  to  be  zero 
[14,31]. 
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(a)  (b) 


Fig.  3.  Microscopic  pictures  of  the  carbon  gas  diffusion  media:  (a)  carbon  cloth  (a  white  square  represents  the  computational  domain  and  a  white  dotted  circle  the 
region  of  tow)  and  (b)  carbon  paper  (a  white  square  represents  the  computational  domain). 


2.4.  LB  model  for  the  Stokes/Brinkman  formulation  in  a 
porous  medium 

Since  the  LB  Stokes/Brinkman  formulation  is  well  explained 
in  Spaid  and  Phelan  [15]  and  Park  et  al.  [16],  here,  the 
formulation  is  briefly  introduced  for  further  extension  to  multi- 
component  LB  model.  Fluid  flow  in  the  open  regions  is  governed 
by  the  Stokes  equation  given  by 

/W2«  =  VP,  (26) 

where  /x  is  viscosity,  u  is  velocity  vector  and  P  is  pressure.  Inside 
the  porous  tows,  the  flow  is  modeled  by  the  Brinkman  equation: 

MV 2  (u)--f--  (u)  =  V  (P) ,  (27) 

*MOW 

where  Ktow  is  the  permeability  tensor  for  the  fiber  bundles.  To 
recover  the  Brinkman  equation  using  a  LB  method,  it  is  nec¬ 
essary  to  modify  the  standard  equilibrium  distribution  function 
to  reduce  the  magnitude  of  the  momentum,  leaving  the  direc¬ 
tion  of  momentum  unchanged.  This  can  be  achieved  by  altering 
the  velocity  u(x,  t)  in  the  equilibrium  distribution  function  by 
incorporating  a  forcing  term  given  by 

zF(x,  t) 

U  =  u(x,t)  +  s(x) — - - — ,  (28) 

P(x,  t) 


where  U  replaces  u  in  Eq.  (5)  for  the  equilibrium  distribution 
function,  and  the  variable  s(x,  t)  is  either  0  or  1  depending  on 
whether  a  lattice  site  is  located  in  a  void  or  porous  region,  respec¬ 
tively.  The  porous  region  indicates  the  fibrous  area  with  flow 
resistance  caused  by  the  friction  on  the  fiber  surface  while  the 
void  region  indicates  the  empty  space  between  fiber  bundles  with 
no  direct  effect  of  friction  as  shown  in  Figs.  3(a)  and  4(a).  The 
form  of  the  forcing  term  F(x,  t)  which  is  needed  to  recover  the 
Brinkman  equation  is  given  by 

F(x,  t)  =  —fip(x,  t)u(x ,  t),  (29) 

where  ft  is  a  parameter  controlling  the  magnitude  of  the  momen¬ 
tum  sink.  From  Eqs.  (28)  and  (29)  the  equilibrium  velocity  in  a 
Brinkman  site  (s(x,  t)  =  1.0)  is  defined  as 

U  =  u(x ,  0(1  -  M.  (30) 

Then,  the  LB  equation  in  a  Brinkman  site  is  given  for  the 
steady  state  as 

9  1 

i N2u-pu  =  -S7P.  (31) 

P 

Eq.  (31)  reproduces  the  Brinkman  equation  if  /3  =  v/Kt0W 
in  which  both  the  kinematic  viscosity  and  the  tow  perme¬ 
ability  are  expressed  in  lattice  units  [20].  The  validity  of  LB 


Fig.  4.  Lattice  Boltzmann  models  for  gas  diffusion  media:  (a)  2D  porous  model  for  carbon  cloth  medium  and  (b)  2D  micro-obstacle  model  for  carbon  paper  (Case 
A). 
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Stokes/Brinkman  modeling  was  verified  by  comparing  the  LB 
simulation  results  with  analytical  and  experimental  measure¬ 
ments  in  a  previous  work  [16]. 

3.  Results  and  discussion 

The  GDL  of  a  PEM  fuel  cell  has  been  numerically  modeled 
using  LB  method  via  two  different  LB  formulations.  The  woven 
carbon  cloth  is  modeled  as  void  space  and  porous  area  which 
has  certain  permeability  by  the  LB  Stokes/Brinkman  equation 
[15,16]  while  randomly  distributed  fibrous  porous  structure  of 
the  carbon  paper  is  modeled  based  on  the  actual  microscopic 
images  using  the  LB  solid  obstacle  modeling.  In  Section  3.1, 
the  permeability  of  the  modeled  medium  is  compared  with 
experimental  measurements  in  literature  [10].  The  effect  of  fiber 
orientation  on  the  permeability  of  the  porous  medium  has  been 
investigated  within  a  typical  range  of  porosity.  In  the  follow¬ 
ing  section,  the  two  porous  LB  flow  models  are  extended  to 
unsteady  multi-phase  LB  simulations  in  the  porous  medium. 
The  LB  model  simulates  unsteady  motion  of  a  liquid  droplet 
and  its  breaking  up  during  its  passage  through  a  complicated 
porous  medium  being  driven  by  ambient  fluid  convection. 

3.1.  LB  models  for  the  gaseous  reactant  flow  in  an 
inhomogeneous  porous  medium 

Fig.  3(a)  and  (b)  presents  the  microscopic  pictures  of  the 
woven  carbon  cloth  and  carbon  paper,  respectively,  which  are 
the  most  common  materials  for  the  gas  diffusion  layers  in  PEM 
fuel  cells.  The  porous  structure  of  the  carbon  cloth  consists  of 
carbon  fiber  tows,  the  bundles  of  carbon  fiber,  and  void  spaces 
among  tows.  The  combinational  effect  of  the  void  space  and  tow 
permeability  results  in  the  effective  permeability  of  the  porous 
medium.  Since  the  fiber  tow  has  much  finer  substructures  inside 
compared  to  its  characteristic  length  scale  (i.e.  the  diameter  of 
the  fiber  tow)  it  would  be  numerically  efficient  to  model  the  fiber 
tow  as  a  uniform  porous  area  with  certain  permeability  neglect¬ 
ing  its  substructure.  In  Fig.  4(a),  the  carbon  cloth  is  modeled  as 
consisting  of  void  space  and  fiber  tow  area  which  has  certain  per¬ 
meability  and  the  LB  Stokes/Brinkman  formulation,  Eqs.  (l)-(5) 
with  extra  forcing  terms  (28)-(30)  in  porous  area,  are  solved  for 
a  transverse  flow  over  2D  porous  circle  [15,16].  The  size  of 
physical  domain  is  400  [xm  x  400  p,m  as  shown  in  Fig.  4(a)  and 
50  x  50  grid  points  are  given  for  the  numerical  domain  to  model 


this  geometry.  The  side  boundaries  are  assumed  to  be  periodic 
based  on  the  fact  that  the  entire  medium  has  the  reciprocating 
pattern  of  the  fiber  tow  and  void  space  in  the  plane  directions. 
Meanwhile,  the  porous  structure  of  the  carbon  paper  may  be 
characterized  as  void  space  surrounded  by  randomly  distributed 
carbon  fibers  with  no  distinct  pattern  found  in  the  carbon  cloth. 
The  Stokes/Brinkman  formulation  is  not  proper  to  model  the 
detail  of  such  complicated  random  structure  which  has  a  dom¬ 
inant  influence  on  the  permeability  of  the  medium.  Therefore, 
in  this  work,  the  true  porous  structure  is  recovered  from  the 
microscopic  image  of  carbon  paper  using  a  large  number  of  grid 
point  (150  x  150)  in  the  physical  domain  of  400  |xm  x  400  p,m 
as  shown  in  Fig.  4(b).  For  this,  the  microscopic  image  of  carbon 
paper  is  converted  into  a  binary  data  (0:  void  and  1:  solid  obsta¬ 
cle)  by  applying  a  certain  value  of  threshold.  The  solid  obstacle 
in  the  numerical  domain  is  an  impermeable  object  with  no  slip 
boundary  condition  on  its  surface  which  is  the  equivalent  sur¬ 
face  characteristic  of  the  carbon  fiber  in  the  physical  domain.  In 
this  way,  the  porosity  of  the  medium  can  be  controlled  accord¬ 
ing  to  the  value  of  the  threshold  without  altering  the  major  fiber 
distribution  and,  in  the  present  LB  solid  obstacle  modeling,  the 
porosity,  0,  is  determined  as  the  ratio  of  the  number  of  void 
grid  points  to  the  number  of  total  grid  points  in  the  converted 
binary  data.  Sample  images  were  taken  at  random  position  of  the 
images  to  investigate  the  effect  of  fiber  orientation  and  porosity 
on  the  resultant  permeability  of  the  medium.  Fig.  5  shows  five 
different  sample  media  of  carbon  paper  which  are  chosen  for  the 
present  LB  simulation.  Medium  A-C  are  characterized  as  dom¬ 
inant  fiber  orientations  with  relatively  low  porosities,  whereas 
medium  D  and  E  are  lack  of  such  dominant  fiber  orientations 
with  relatively  higher  porosities. 

Fig.  6(a)  presents  the  velocity  and  pressure  distributions  in 
the  porous  medium  modeled  by  Stokes/Brinkman  formulation 
when  the  tow  permeability,  K^,  is  10“ 10  m2.  Constant  pressure 
is  applied  at  the  bottom  and  top  surfaces  to  drive  the  flow  motion 
while  the  two  vertical  surfaces  are  assumed  to  be  periodic.  The 
velocity  vectors  are  normalized  by  the  maximum  velocity  in  the 
computational  domain  where  the  tow  area  is  marked  as  a  dashed 
circle.  It  is  seen  that  the  magnitude  of  velocity  is  substantially 
smaller  inside  the  tow  area  than  outside  the  tow  area  due  to 
the  flow  resistance  caused  by  the  low  permeability  of  the  tow. 
The  pressure  contours  are  symmetrical  profiles  for  the  upstream 
and  downstream  divided  by  a  straight  contour  that  is  normal 
to  the  mean  flow  direction  (y-direction).  Fig.  6(b)  shows  veloc- 


Fig.  5.  Numerical  models  of  gas  diffusion  layer  obtained  from  the  various  samples  of  microscopic  images  (carbon  paper);  f  is  the  porosity,  determined  as  the  ratio 
of  the  number  of  void  grid  points  to  the  number  of  total  grid  points  in  the  converted  binary  data. 
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Fig.  6.  Lattice  Boltzmann  simulation  results  in  the  gas  diffusion  media;  the  vector  represents  the  fluid  velocity  and  the  color  represents  the  pressure  distribution: 
(a)  porous  model  for  carbon  cloth  medium,  ^Ttow  =  10“ 10  m2,  (b)  micro-obstacle  model  for  carbon  paper;  black  colored  area  represents  solid  obstacle  (Model  A), 
Rec  -  0.05  and  the  pressure  and  velocity  vectors  are  non-dimensionalized  using  average  pressure  and  maximum  velocity  in  the  domain. 


ity  and  pressure  distributions  in  a  porous  medium  modeled  by 
the  LB  solid  obstacle  formulation  for  the  Medium  A  in  Fig.  5. 
The  black  colored  area  represents  the  impermeable  solid  obsta¬ 
cles  (carbon  fibers)  and  it  shows  a  complicated  flow  field  in  the 
void  space  among  fibers.  The  flow  is  squeezed  at  several  narrow 
paths  resulting  in  a  high  velocity  with  a  rapid  pressure  drop. 
It  is  also  observed  that  a  major  flow  blockage  can  be  caused  by 
small  obstacles  regularly  distributed  in  a  finite  area  as  seen  in  the 
left  upper  part  of  Fig.  6(b).  On  the  contrary,  the  large  obstacles 
parallel  with  the  main  flow  direction  do  not  necessarily  pro¬ 
duce  large  flow  resistance.  This  implies  that  the  permeability  of 
the  porous  medium  is  not  only  dependent  on  the  porosity  but 
also  have  strong  influence  from  the  characteristics  of  the  porous 
structure. 

To  investigate  the  effect  of  the  porous  structure  such  as  fiber 
orientations  on  the  permeability  of  the  porous  medium,  the  flow 
simulation  is  performed  for  the  two  major  perpendicular  direc¬ 
tions  for  the  various  samples  of  porous  medium  in  Fig.  5.  The 
permeability  of  the  medium  is  calculated  based  on  the  Darcy’s 
law: 


M«) 

V(P>’ 


(32) 


where  (u)  and  ( P )  are  the  average  velocity  and  pressure  in  the 
computational  domain,  respectively.  The  value  of  permeability 
is  non-dimensionalized  by  the  length  of  computational  domain, 
L.  In  Fig.  7,  Medium  A  and  C  have  higher  permeability  when 
the  pressure  gradient  is  given  for  the  y-direction  due  to  the  fact 
that  main  fiber  orientations  is  also  parallel  with  the  y-direction  as 
shown  in  Fig.  5.  Similarly,  the  Medium  B  has  higher  permeabil¬ 
ity  for  the  v-direction  with  the  main  fiber  direction  parallel  with 
v-direction.  However  such  phenomenon  is  not  apparent  in  case 
of  Medium  D  and  E  since  the  fiber  orientation  is  not  obvious  for 
those  two  samples.  The  value  of  permeability  varies  in  a  similar 
range  with  that  of  Medium  A-C  except  for  the  cases  of  A-x  and 
C-x  in  Fig.  7  although  the  porosity  of  Medium  D  and  E  are  quite 


higher  than  other  three  media.  It  seems  that  the  effect  of  the  fiber 
orientation  on  the  permeability  of  the  medium  can  be  dominant 
over  that  of  porosity  within  the  given  range  of  porosity. 

Fig.  8  compares  the  permeability  of  the  porous  medium 
obtained  by  the  present  LB  solid  obstacle  modeling  with  the 
experimental  measurements  of  in-plane  permeability  for  the  two 
different  GDLs  in  Feser  et  al.  [10] .  SGL  3 1BA  (SGL  Carbon)  and 
TGP-60-H  (Toray)  are  a  non-woven  carbon  fiber  and  a  paper- 
based  carbon  fiber,  respectively,  without  micro-porous  layers. 
Both  samples  contain  almost  identical  porous  structures  con¬ 
sisting  of  randomly  oriented  carbon  fibers  as  the  one  shown  in 
Fig.  5(b).  The  porosity  of  Medium  A  is  modified  by  adjusting  the 
value  of  threshold  during  the  image  converting  process.  The  LB 
simulation  results  present  two  discrete  lines  for  the  different  flow 
directions  in  Fig.  8  however  each  line  does  not  show  much  incre¬ 
ment  as  the  porosity  of  the  medium  increases  from  70  to  80%. 
This  is  attributed  to  the  fact  that  the  fibers  become  thinner  but  the 
orientations  are  not  altered  as  the  porosity  is  increased.  It  is  also 
observed  that  the  experimental  data  for  various  gas  diffusion  lay- 
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Fig.  7.  Permeability  of  the  porous  medium  for  the  two  main  flow  directions:  x 
and  y  in  the  legend  indicate  the  flow  direction  considered  in  the  simulation. 
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Fig.  8.  Comparison  of  the  permeability  of  the  porous  medium  (Medium  A)  with 
experimental  measurements  of  in-plane  permeability  of  various  gas  diffusion 
layers  (carbon  papers)  in  Feser  et  al.  [10]. 


ers  in  Feser  et  al.  [10]  fall  between  these  two  principal  values  of 
permeability.  Although  it  is  not  so  meaningful  to  compare  the  LB 
Stokes/Brinkman  modeling  and  the  experimental  measurements 
since  the  correlation  between  the  porosity  and  permeability  of 
the  fiber  tow  is  not  clearly  known  for  the  present  porous  medium, 
the  permeability  of  the  porous  medium  can  be  calculated  as 
2.05  x  10_10m2  using  LB  Stokes/Brinkman  model  when  the 
diameter,  porosity  and  permeability  of  the  fiber  tow  are  assumed 
to  be  360  pum,  0.65  and  1.23  x  10“ 10  m2,  respectively,  so  that  the 
medium  has  the  porosity  of  0.77.  This  value  is  one  order  of  mag¬ 
nitude  larger  than  the  experimental  measurements  in  Fig.  8  and 
the  difference  is  mainly  attributed  to  neglecting  the  structural 
effect  in  the  Stokes/Brinkman  model. 


3.2.  LB  multi-phase  flow  simulation  in  a  porous  medium 

One  of  the  primary  advantages  of  LB  methods  for  simulating 
fluid  flows  as  compared  to  the  traditional  numerical  methods  is 
their  ability  to  robustly  model  interfaces  between  two  or  more 
fluids  in  complicated  geometries  such  as  porous  medium.  The 
mechanism  of  liquid  water  accumulation  in  the  GDL  is  of  great 
interest  for  PEM  fuel  cell  researchers  since  the  phenomenon  is 
widely  accepted  as  one  of  the  main  mechanisms  which  limit  the 
performance  of  a  PEM  fuel  cell.  In  this  section,  the  two  flow 
models  in  the  previous  section  are  extended  to  multi-phase  LB 
simulations  by  incorporating  the  inter-particle  interaction  model 
[20,29].  The  magnitude  of  inter-particle  force  is  set  through  the 
value  of  gaa  in  Eq.  (10).  Then  the  viscosity  of  the  fluid  is  given 
as 


where  /3C r  is  the  mass  density  concentration  of  the  kth  component 
and  is  defined  as 

Pa  =  (34) 

Z^crPo- 

The  multi-phase  flow  in  porous  medium  involves  low  velocity 
and  dominant  effect  of  capillary  forces  between  fluid/fluid  and 
fluid/surface  interfaces.  The  phenomenon  can  be  characterized 
by  the  capillary  number  which  represents  the  relative  effect  of 
viscous  forces  versus  surface  tension  acting  across  an  interface 
between  a  liquid  and  a  gas,  or  between  two  immiscible  liquids. 


Fig.  9.  Unsteady  two-phase  simulation  results  using  the  LB  Stokes-Brinkman  model  incorporated  with  inter-particle  interaction  model  [29];  blue  area  is  the  liquid 
water  and  the  interior  of  the  dotted  circle  indicates  the  porous  region  representing  a  fiber  tow  of  the  carbon  cloth  GDL.  (1  and  2)  Liquid  droplet  starts  to  move,  driven 
by  the  mean  flow  motion;  (3-6)  droplet  deforms  and  stretched  around  the  porous  region;  (7  and  8)  detached  droplets  merge  and  leave  the  porous  area  while  a  part 
of  liquid  water  is  trapped  and  remains  inside  in  the  porous  area.  Re l  =  0. 1 ,  Ca  =  2.3  x  10-5  and  ^Ttow  =  10-9  m2 .  (For  interpretation  of  the  references  to  color  in  this 
artwork,  the  reader  is  referred  to  the  web  version  of  the  article.) 
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Fig.  10.  Unsteady  two-phase  simulation  results  using  LB  binary  fluid  model  simulating  two-phase  flow  through  fibrous  structure  of  carbon  paper  GDL;  blue  area 
is  the  liquid  drop  and  black  colored  area  indicates  solid  obstacle  (carbon  fiber).  (1  and  2)  Liquid  droplet  starts  to  move,  driven  by  the  mean  flow  and  hit  the  porous 
surface;  (3-8)  droplet  deforms  and  breaks  up  into  smaller  ligaments  moving  into  the  porous  region;  (9-12)  detached  droplets  leave  the  porous  area  while  the  major 
part  of  liquid  water  is  trapped  and  remains  in  the  porous  area.  Re^  =  0.1  and  Ca  =  2.3  x  10-5.  (For  interpretation  of  the  references  to  color  in  this  artwork,  the  reader 
is  referred  to  the  web  version  of  the  article.) 


It  is  defined  as 


£ 


(35) 


where  V  is  a  characteristic  velocity  and  £  is  the  surface  or  inter¬ 
facial  tension  between  the  two  fluid  phases.  In  this  work,  the 
surface  tension  between  each  fluid  is  estimated  from  Laplace’s 
law  following  Kang  et  al.  [34]  and  the  resultant  capillary  num¬ 
ber,  Ca ,  for  this  simulation  is  found  to  be  2.3  x  10-5  that  is 
within  a  typical  range  of  Ca  for  the  air/liquid  water  interface  in 
a  porous  medium.  In  case  of  Stokes/Brinkman  model,  the  veloc¬ 
ity  of  each  component  is  replaced  with  Eq.  (28)  in  the  porous 
region  as  a  result  of  drag  force.  The  Reynolds  number  is  defined 
using  the  average  velocity  at  the  inlet,  V,  and  the  thickness  of 
the  computational  domain,  L  as 

ReL  =  —.  (36) 

I1 

Present  model  consists  of  two  fluids  (blue  and  white), 
as  shown  in  Fig.  9,  which  may  represent  liquid  water  and 
gaseous  reactant  in  a  PEM  fuel  cell,  respectively.  Fig.  9  presents 
unsteady  multi-component  LB  simulation  results  based  on  the 


Stokes/Brinkman  formulation.  The  relaxation  time,  r,  is  fixed 
as  unity  for  both  components  and  the  density  ratio  of  the  blue 
to  white  fluid  is  set  to  be  2.0.  The  top  and  bottom  surfaces  are 
assigned  to  be  no  slip  walls  and  the  constant  pressure  is  given 
for  the  white  fluid  at  the  inlet  (left  wall)  and  outlet  (right  wall) 
boundaries  to  push  the  blue  fluid  toward  the  v-direction.  As  a 
result,  the  blue  droplet  starts  to  move,  driven  by  the  mean  flow 
motion  of  the  white  fluid  as  shown  in  Fig.  9(1)  and  (2).  The  blue 
droplet  deforms  continually  and  is  stretched  around  the  porous 
region  pushed  by  the  white  fluid  in  Fig.  9(3)-(5).  In  the  mean¬ 
while,  the  stretched  part  of  blue  fluid  penetrates  into  the  porous 
region.  The  blue  fluid  eventually  breaks  up  into  three  subparts 
in  Fig.  9(6).  The  two  detached  droplets  flow  around  the  porous 
region  and  then  merge  together  as  one  larger  droplet  in  Fig.  9(7). 
This  larger  droplet  continues  to  move  along  with  the  mean  flow 
motion  while  the  smaller  one  left  behind  is  trapped  in  the  porous 
area  as  shown  in  Fig.  9(8). 

Fig.  10  presents  another  unsteady  two-phase  simulation 
of  the  liquid  droplet  passing  through  the  complicated  porous 
medium  representing  carbon  paper  GDL  (Medium  A)  using  LB 
solid  obstacle  modeling  incorporated  with  multi-component  LB 
model.  The  boundary  conditions  are  the  same  as  those  of  previ- 
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ous  two-phase  LB  simulation.  The  surface  condition  of  fibrous 
porous  structure  is  controlled  by  the  value  of  Gaw  in  Eq.  (12), 
and  is  given  as  0.1  to  generate  wetting  condition  for  the  present 
case.  It  should  be  noted  that  the  mixed  wetting  surface  condition 
[35-37]  in  a  practical  GDL  may  be  modeled  by  alternating  the 
value  of  Gaw  for  selected  surfaces  or  regions  in  the  domain.  In 
Fig.  10(1)  and  (2),  it  is  seen  that  the  large  blue  droplet  driven 
by  the  mean  flow  motion  hits  the  porous  surface.  As  the  droplet 
is  pushed  by  the  white  fluid  it  starts  to  deform  and  is  squeezed 
into  the  void  space  in  the  porous  medium  in  Fig.  10(3)  and 
(4)  and,  later  a  chunk  of  the  blue  fluid  is  detached  and  accu¬ 
mulated  in  the  upper  part  of  porous  medium  in  Fig.  10(5). 
The  major  part  of  blue  fluid  breaks  up  into  three  subparts  in 
Fig.  10(6)-(8)  and  each  detached  sub-droplets  reshape  its  form 
due  to  the  surface  tension  of  the  fluid.  In  Fig.  10(9)-(12),  the 
detached  droplets  leave  the  porous  area  while  other  three  sub¬ 
parts  of  liquid  water  are  trapped  and  remain  in  the  void  space  of 
the  medium.  Present  LB  model  predicts  well  the  unsteady  liq¬ 
uid  water  accumulation/removal  process  under  a  strong  ambient 
convective  flow  in  a  porous  medium.  The  model  would  be  useful 
tool  to  investigate  the  liquid  water  removal  from  a  gas  diffusion 
layer  under  land  area  utilizing  cross-leakage  flow  [9] .  Although 
not  shown  here,  higher  pressure  and  hydrophobic  surface  con¬ 
dition  resulted  in  more  efficient  liquid  water  removal  from  the 
porous  medium.  While,  for  strict  comparison  with  experimental 
observations,  it  may  require  higher  density  ratio  between  each 
fluid  which  is  the  case  of  the  interface  between  air  and  liquid 
water.  The  density  ratio  of  the  present  multi-component  model 
is  confined  to  0(1)  due  to  the  stability  issues  [20,29].  There 
exist  a  few  recently  developed  LB  models  [38,39]  available  for 
the  multi-phase  flow  simulation  in  a  higher  density  ratio  up  to 
1000  although  their  stability  is  not  yet  proven  in  a  complicated 
geometry. 

4.  Concluding  remarks 

In  this  work,  lattice  Boltzmann  method  has  been  applied  to 
model  the  multi-phase  micro-scale  flow  through  the  gas  diffu¬ 
sion  layers  in  a  PEM  fuel  cell  that  are  made  of  carbon  cloth 
and  carbon  paper,  respectively.  The  woven  carbon  cloth  is  mod¬ 
eled  as  void  space  and  porous  area  with  certain  permeability 
while  randomly  distributed  fibrous  structure  of  the  carbon  paper 
is  directly  modeled  based  on  the  microscopic  images  of  the 
structure  and  it  is  treated  as  impermeable  obstacles.  The  model 
predictions  match  well  with  the  experimental  measurements  of 
permeability  in  literature.  The  effect  of  fiber  orientation  on  the 
permeability  of  the  medium  is  investigated  on  the  basis  of  var¬ 
ious  sample  images  of  carbon  paper.  The  permeability  varies 
considerably  according  to  the  flow  direction  indicating  that  the 
effect  of  fiber  orientation  is  critical  to  the  directional  dependence 
of  the  permeability  of  inhomogeneous  porous  medium.  Multi¬ 
phase  lattice  Boltzmann  porous  model  has  been  developed  on 


the  basis  of  the  two  LB  flow  models  by  incorporating  multi- 
component  LB  model  with  inter-particle  forces.  The  model 
successfully  simulates  the  unsteady  characteristic  behaviors  of 
liquid  droplet  motion  in  the  porous  medium  providing  a  useful 
tool  to  investigate  the  mechanism  of  liquid  water  accumula¬ 
tion/removal  in  the  porous  electrode  of  a  PEM  fuel  cell. 

Acknowledgment 

This  work  was  supported  by  AUT021,  the  Network  of  Cen¬ 
ters  of  Excellence,  Canada. 

References 

[1]  V.A.  Paganin,  E.A.  Ticcianelli,  E.R.  Gonzalez,  J.  Appl.  Electrochem.  26 
(3)  (1996)  297. 

[2]  J.  Benziger,  J.  Nehlsen,  D.  Blackwell,  T.  Brennan,  J.  Itescu,  J.  Membr.  Sci. 
261  (2005)  98-106. 

[3]  G.  Inoue,  Y.  Matsukuma,  M.  Minemoto,  J.  Power  Sources  154  (2006)  8-17. 

[4]  H.  Lee,  J.  Park,  D.  Kim,  T.  Lee,  J.  Power  Sources  131  (2004)  200-206. 

[5]  L.  Jordan,  A.  Shukla,  T.  Beehrsing,  N.  Avery,  B.  Muddle,  M.  Forsyth,  J. 
Power  Sources  (2000)  250-254. 

[6]  J.  Park,  X.  Li,  J.  Power  Sources  163  (2007)  853-863. 

[7]  G.  Inoue,  Y.  Matsukuma,  M.  Minemoto,  J.  Power  Sources  157  (2006) 
136-152. 

[8]  T.  Kanezaki,  X.  Li,  JJ.  Baschuk,  J.  Power  Sources  162  (2006)  415-425. 

[9]  J.  Park,  X.  Li,  D.  Tran,  T.  Abdel-Baset,  D.S.  Hussey,  D.L.  Jacobson,  M. 
Arif,  Int.  J.  Hydrogen  Energ.,  in  press. 

[10]  J.P.  Feser,  A.K.  Prasad,  S.G.  Advani,  J.  Power  Sources  162  (2006)  1226. 

[11]  S.  Dutta,  S.  Shimpalee,  J.W.  Van  Zee,  Int.  J.  Heat  Mass  Transf.  44  (2001) 
2029. 

[12]  J.G.  Pharoah,  J.  Power  Sources  144  (2005)  77. 

[13]  J.  Park,  X.  Li,  J.  Power  Sources  163  (2007)  853. 

[14]  A.S.  Joshi,  K.N.  Gre,  A.A.  Peracchio,  W.K.S.  Chiu,  J.  Power  Sources  164 
(2007)  631. 

[15]  M.  Spaid,  F.R.  Phelan  Jr.,  Phys.  Fluids  9  (1997)  2468. 

[16]  J.  Park,  M.  Matsubara,  X.  Li,  J.  Power  Sources  173  (2007)  404-414. 

[17]  F.R.  Phelan  Jr.,  G.  Wise,  Composites  27A  (1995)  25. 

[18]  X.  He,  L.  Luo,  J.  Stat.  Phys.  88  (1997)  927. 

[19]  A.  Cali,  S.  Succi,  A.  Cancelliere,  R.  Benzi,  M.  Gramignani,  Phys.  Rev.  A 
45  (1992)  5771. 

[20]  X.W.  Shan,  H.D.  Chen,  Phys.  Rev.  E  47  (1993)  1815. 

[21]  B.J.  Palmer,  D.R.  Rector,  Phys.  Rev.  E  61  (2000)  5295. 

[22]  S.P.  Dawson,  S.  Chen,  G.D.  Doolen,  J.  Chem.  Phys.  98  (1993)  1514. 

[23]  X.  He,  N.  Li,  Comput.  Phys.  Commun.  129  (2000)  158. 

[24]  B.  Li,  D.Y.  Kwok,  Int.  Heat  Mass  Transf.  46  (2003)  4235. 

[25]  J.  Park,  K.Y.  Huh,  X.  Li,  J.  Electroanal.  Chem.  591  (2006)  141. 

[26]  PL.  Bhatnagar,  E.P.  Gross,  M.  Krook,  Phys.  Rev.  94  (1954)  511. 

[27]  X.  He,  L.  Luo,  Phys.  Rev.  E  56  (1997)  6811. 

[28]  H.  Chen,  S.  Chen,  W.H.  Matthaeus,  Phys.  Rev.  A  45  (1992)  R5339. 

[29]  X.  Shan,  H.  Chen,  Phys.  Rev.  E  49  (1994)  2941. 

[31]  Q.  Zou,  X.  He,  Phys.  Fluids  9  (1997)  1591-1598. 

[34]  Q.  Kang,  D.  Zhang,  S.  Chen,  Phys.  Fluids  14  (2002)  3203. 

[35]  E.C.  Kumbur,  K.V.  Sharp,  M.M.  Mench,  J.  Power  Sources  168  (2007)  356. 

[36]  N.  Djilali,  D.  Lu,  Int.  J.  Therm.  Sci.  41  (2002)  29. 

[37]  Q.  Ye,  T.V.  Nguyen,  J.  Electrochem.  Soc.  154  (2007)  B1242. 

[38]  T.  Inamuro,  T.  Ogata,  S.  Tajima,  N.  Konishi,  J.  Comput.  Phys.  198  (2004) 
628. 

[39]  T.  Lee,  C.  Lin,  J.  Comput.  Phys.  206  (2005)  16. 


